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Abstract 



We describe a variational method to solve the Holstein model for an elec- 
tron coupled to dynamical, quantum phonons on an infinite lattice. The 
variational space can be systematically expanded to achieve high accuracy 
with modest computational resources (12-digit accuracy for the Id polaron 
energy at intermediate coupling). We compute ground and low- lying excited 
state properties of the model at continuous values of the wavevector k in 
essentially all parameter regimes. Our results for the polaron energy band, 
effective mass and correlation functions compare favorably with those of other 
numerical techniques including DMRG, Global Local and exact diagonaliza- 
tion. We find a phase transition for the first excited state between a bound 
and unbound system of a polaron and an additional phonon excitation. The 

phase transition is also treated in strong coupling perturbation theory. 
PACS: 71.38.+i, 72.10.Di, 71.35.Aa 
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I. INTRODUCTION 



Polaron formation as a consequence of electron-phonon coupling appears in many con- 
texts in condensed matter physics, including charge carriers in colossal magnetoresistance 
materials and in high-temperature superconductors As theoretical research in this 
field spans over five decades, many analytical techniques have been applied to this problem 
||. The applicability of these methods is usually limited to a particular parameter regime, 
frequently far from the physically most interesting crossover regime. Despite extensive an- 
alytical work in this field, there remain many open problems including the nature of the 
crossover to large polaron mass, the form of various correlation functions, the nature of the 
polaron excited states, and dynamic properties of the polaron. 

Constantly growing computer capabilities have allowed research in this field to take 
advantage of various numerical methods, including exact diagonalization techniques (ED) 
H 12] , quantum Monte Carlo calculations (QMC) fl3| , |14]| , variational methods including the 
Global-Local method (GL) [T5|,[16j, and recently developed density matrix renormalization 
group techniques (DMRG) fl7|| . A comparison of results obtained by different methods for 



energy bands and effective masses is contained in the work of Romero et al. [15|. Although 
most of these methods give reliable results in a wide range of parameter regimes, each suffers 
from different shortcomings. The ED technique gives reliable results on small lattices (up 
to 20 sites J8|) for ground and excited state properties. Results are limited to discrete mo- 
mentum A;— points. QMC methods can treat large system sizes (over 1000 lattice sites) and 
provide accurate results for thermodynamic polaron properties. Dynamic properties, how- 
ever, require analytic continuation from imaginary time, which is an ill-posed problem that 
is extremely sensitive to statistical noise. The GL method gives reliable results for energy 
bands and the polaron effective mass on reasonably large system sizes (32 sites), however, 
it is limited to ground state properties. The DMRG method seems to be most successful in 
computing ground state properties including various correlation functions. Despite the lack 
of translational symmetry, it provides reasonably accurate results for energy bands and the 
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effective mass. It can deal with large system sizes (e.g. 80 sites and 30 phonons per site), 
delivering reliable results in the intermediate and strong coupling regime. DMRG is less 
reliable in the weak coupling regime. It generally treats large systems with open boundary 
conditions and does not allow the calculation of dynamic or large k properties or excited 
states, although there are some exceptions to these rules [|T^|]. 

In this paper we focus on several issues of the Holstein problem. We present a simple and 
computationally efficient method based on the exact diagonalization technique. We apply 
the method to the one- dimensional, single-electron Holstein model. Among the advantages 
of our method are the simplicity of our approach, the efficiency in selecting the variational 
space, and the ability to compute both ground and excited state physical properties of the 
system at continuous total wavevector k. We define the variational space on an infinite 
lattice. Even though most of our calculations were done on a workstation, our results are 
often superior to those of other numerically more intensive methods fl8|,|9|,[T7|j . We test the 
method by comparing results for the energy band, effective mass, quasiparticle weight and 
correlation functions with some of the most successful recent numerical methods. 

The second part of this paper is devoted to the investigation of the first excited state of 
the model. Using our numerical method as well as a strong coupling analytical approach, 
we find a phase transition between a state where the polaron forms a bound state with an 
additional phonon excitation, and one where there is no bound state. 

Generalizations of this method can even be used to calculate coherent quantum dynamics 
far from equilibrium, including those of a polaron driven by a strong electric field fl9| , and 



of an electron tunneling across a potential drop while coupled to phonons |20 
We consider the Holstein Hamiltonian JIT] 



H=-t^2(4c j+1 + H.c.) 

j 
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where cj creates an electron and a] creates a phonon on site j. We consider the case where 
a single electron hops between nearest neighbor Wannier orbitals in Id, and interacts with 



dispersionless optical phonons. The electron-phonon coupling strength is A, and is local in 
real space (/c-independent). The parameters t, u, and A all have units of energy, and can 
be used to form two dimensionless ratios. (Different conventions for the parameters are 
sometimes used in other papers.) 

It is clear at the outset that for any finite values of the parameters, the exact groundstate 
will be a delocalized state with momentum k, and not a "self-trapped" solution that breaks 
translation-invariance. The simple argument is that if the ground state were self-trapped, its 
center could be shifted and a degenerate ground state obtained. If the Hamiltonian has any 
nonzero matrix element between states with different centers, a lower energy state can be 
obtained by a phased superposition of wavefunctions with different centers (a Bloch state). 
The only known way that this argument can fail is with (unphysically) strong electron- 
phonon coupling as uo — > to a gapless phonon spectrum, in which case the matrix element 
between different centers can vanish |22|] . In that case, an exact ground state can be written 
as either self-trapped or as a Bloch state. Self-trapping cannot occur here because the 
phonon spectrum has a gap. Although the exact ground state is a Bloch state, for some 
parameters (e.g. very large A), a self-trapped solution can have a low variational energy [23 



It has indeed been proven that polaron ground state properties, including the energy and 



effective mass, are analytic functions of the Hamiltonian parameters |24| . We will see below 
that such a theorem cannot be extended to include excited states. 



II. THE METHOD 

A complete set of basis states for the many-body Hilbert space can be written 

\M) = \j;...,n j - 1 ,n j ,n j+ i,...), (2) 

where j is the electron site and there are n m phonons on site m. A variational subspace is 
constructed beginning with an initial state, taken to be an electron on site j = with no 
phonons, and operating repeatedly (iVVtimes) with the off-diagonal pieces of the Hamilto- 
nian, Eq. (p. At each step, a basis state is added when there is a nonzero t or A matrix 
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element to a state previously in the space. These states and all of their translations on an 
infinite lattice are included in the variational space. (A translation moves the electron and 
all phonon excitations m sites to the right.) If a basis state can be generated in more than 
one way, only one copy is retained. All nonzero matrix elements of the Hamiltonian between 
retained basis states are included. 

A small variational Hilbert space is shown in Fig. (|l|). The dots may be thought of 
as basis states in the many-body Hilbert space, or alternatively as Wannier orbitals in a 
periodic (one-body) tight-binding model, with hopping matrix elements given by the bonds. 
The variational space shown in Fig. (|l|) is still infinite, since the lattice repeats periodically 
to infinity. It is clear from Bloch's theorem, however, that each eigenstate can be written 
as e J ' 1 S r m , where k is the momentum, j is the unit cell, and fy m is a set of M complex 
amplitudes, one for each state in the unit cell (M = 7 for the example of Fig. [I]). For a 
given momentum k, the resulting numerical problem is to find the eigenstates of an M x M 
(sparse) hermitian matrix using Lanczos or another method. 

Figure (|2[) plots the energy eigenvalues for a larger variational space containing a maxi- 
mum of 9 phonon excitations. The figure superficially resembles a "band structure" , which 
however encodes ground and excited state information for the many-body (many phonon) 
polaron problem. The AC conductivity of the polaron, for example, appears as an "inter- 
band" transition in this mapping. 

The largest variational basis that we have used has Nh = 22, or M = 1.2 x 10 7 states. 
It is usually unnecessary to use such a large basis for intermediate-coupling ground state 
properties, even to obtain 12-digit accuracy for the energy. The variational Hilbert space 
we construct is not a standard one, and appears to add basis states more efficiently than 
some other methods. A basis state is included if it can be reached using N\ phonon creation 
operators and N t electron hops in any order with N t + N\ < N^. For a given Nh, there 
is a basis state with Nh phonon quanta on the same site as the electron and no phonon 
excitations elsewhere. There is also a basis state with Nh — 1 quanta on the site adjacent 
to the electron and no phonon excitations elsewhere. For Nh odd, there is a basis state 
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with (Nh — l)/2 quanta on the electron site and an equal number simultaneously on a 
first neighbor site. The maximum distance that a phonon excitation can appear from the 
electron is I = Nh — 1, but then only a single quantum and only if there are no phonons 



excited elsewhere in the system ||25|| . Only the basis states from a single unit cell (a single 
electron position) are stored in computer memory. 

The energy of a polaron for a finite chain of N sites with periodic boundary conditions 
is lower than that for an infinite lattice with the same parameters. (This is easily to verify 
in weak-coupling perturbation theory, or numerically.) Thus, previous exact diagonalization 
and most other variational approaches produce energies that are variational for the particular 
lattice size N that they treat, but are not variational in the thermodynamic limit iV — > oo. 
The energy we calculate is, in contrast, variational in the thermodynamic limit. (The quoted 
DMRG energies are extrapolated and not variational, although they are generally rather 
accurate.) 

Having the capability to compute the polaron energy E(k) at any k rather than being 
limited to multiples of 2tt/N makes our method more accurate for computing the effective 
mass of the polaron using the standard formula 

m 1 d 2 E(k) 



fc=o> 



(3) 



m* 2t dk 2 

where mo = 1/ (2t) is the effective mass of a free electron. The second derivative is evaluated 
by small finite differences in the neighborhood of k = 0. Note that although the calculated 
energy E(k) is a variational bound for the exact energy, there is no such control on the mass, 
which may be either above or below the exact answer, and is expected to be more difficult 
to obtain accurately. Nevertheless, in the intermediate coupling regime where our method 
at Nh = 20 gives an energy accuracy of 12 decimal places, one can calculate the effective 
mass extremely accurately (6-8 decimal places) by letting Ak — ► 0. 

Further information about the quasiparticle may be obtained by computing the quasi- 
particle residue, the overlap (squared) between a bare electron and a polaron, 

^ = |(V> fe |4lo>l 2 , (4) 



where |0) is the state with no electron and no phonon excitations, and \ipk) is the polaron 
wavefunction at momentum k. The numerical results for Zk =0 given in the next section 
differ by less than 1% from results for the reciprocal effective mass mo/m* obtained from 
Eq. (^). At finite k, Zk provides information about the electronic character of the polaronic 
state. The phonon contribution to the quasiparticle can be measured by the k— dependent 
mean phonon number 

< = EK^l4a# fc )| 2 - (5) 

i 

To describe the polaron at different k we have also computed the static correlation 
function between the electron position and the oscillator displacement 

x (i-j) = {i) k \c\ci((ij + at)|V>fc), (6) 

and the distribution of the number of excited phonons in the vicinity of the electron 

7(* - j) = (ipklclctalajli/jk). (7) 



III. GROUND STATE RESULTS 

In this section, we compare our results for the ground state properties of the Holstein 
model to those obtained by other numerical methods. We also calculate polaron correlation 
functions at finite values of the Bloch wavevector k. 

We start by comparing energy bands E(k) for two different sets of parameters. Figure 
(0) compares the present method to the DMRG [ I5|,17l , GL |1| and the finite cluster ED 



technique ||. Using = 20 we achieve an accuracy in the thermodynamic limit of 12 
decimal places for small k and at least 4 decimal places for large k. Our results for E(k) 
are presented as continuous curves. Agreement of our results with DMRG and ED is good, 
however there is a slight disagreement with the Global-Local method at larger k, which also 
disagrees with the DMRG method. There is also a slight disagreement with the ED at larger 
k which we attribute to the smaller system size used in finite cluster ED calculations. As 



we will demonstrate later in this paper, the extent of the lattice deformation (the size of 
the polaron) increases as k approaches the Brillouin-zone boundary, which makes finite-size 
calculations more susceptible to finite size effects. 

In table (|) we compare the polaron ground state energy at k = for two different 
parameter sets obtained by our and three other numerical methods, exact diagonalization 
(ED) |p6|| , DMRG [|17[], and Global-Local variational |L5]. Comparisons with a greater 



variety of methods can be found in Romero et al. |15| . We have limited our comparison 
to those methods that are accurate to at least 4 decimal places. Our method converges 
to all of the digits shown using Nh = 15 or M = 88052 basis states for A = u) — 1, and 
N h = 18 or M = 731027 basis states for A = y/2, uj — 1. The ground state energy appears to 
converge exponentially with N h , with the accuracy improving by approximately one order 
of magnitude as Nh is increased by 1. The larger (Nh = 18) calculation runs in under a 
minute on a modest workstation. 

Figure (|J) shows our results for the effective mass Eq. (|3]) computed with Nh = 20 in 
comparison with GL and DMRG methods. The parameters span different physical regimes 
including weak and strong coupling (respectively small and large A/a;), and adiabatic (u/t -C 
1) and antiadiabatic (ui/t ^> 1) regimes. We find good agreement with GL away from strong 
coupling and good agreement in all regimes with DMRG. DMRG calculations are not based 
on finite-fc calculations due to a lack of periodic boundary conditions, so they extrapolate the 
effective mass from the ground state data using chains of different sizes, which leads to larger 
error bars and more computational effort. Notice that their discrete data is slightly scattered 
around our curves. Nevertheless, both methods agree well. A closer look at effective masses 
in the antiadiabatic regime (uj/t = 5.0) and X/u < 2 reveals a slight systematic disagreement 
of our results with both methods in comparison. We have compared our results for effective 
mass obtained on different systems from Nh = 16 with M = 178617 states to Nh = 20 with 
M = 2975104 states and obtained convergence of results to at least 4 decimal places in all 
parameter regimes presented in Fig. (^). Our error is therefore well below the linewidth. 
Even though there is no phase transition in the ground state of the model, the polaron 
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becomes extremely heavy in the strong coupling regime. The crossover to a regime of large 
polaron mass is more rapid in adiabatic regime (smaller uj/t). 

Figure (||) shows the quasiparticle residue Zk and the mean phonon number N% h as 
a function of k for the case of small Ao = X 2 /2ujt (A 2 = 0.4, uj = 0.8) and large Ao 
(A 2 = 3.2, cu = 0.8). The two sets of parameters correspond to the large and small po- 
laron regime respectively ||. The DMRG cannot straightforwardly compute this quantity, 
and we compare our results with cluster calculations. Open symbols represent the results of 
Wellein and Fehske |§ obtained on a N = 14 site cluster for the same choice of parameters. 
Except for the fact that their results are limited to discrete k— points defined by the size of 
their system, we find excellent agreement between the two methods in the weak coupling 
case. In the strong coupling regime there is an approximately 1% disagreement in N% due 
to a lack of phonon degrees of freedom in the variational space of the ED calculation. Our 
results in the weak coupling case show a smooth crossover from predominantly electronic 
character of the wavefunction for small k (large and small Njf 1 ~ 0) to predominantly 
phonon character around k = ir characterized by Z^ pa and N][ pa 1. In the strong 
coupling regime there is less k— dependence. The Z^ is already close to zero at small k, 
indicating strong electron-phonon interactions that lead to a large polaron mass. 

Jeckelmann and White have calculated the electron-phonon correlations at k — for the 
Id and 2d polaron using the DMRG method |i~7| . We compare our results with theirs for the 
Id case. There is good agreement in Fig. (|6]a) for intermediate coupling parameters to = 1 
and A = 0.5. Figure (||d) plots correlations for uj = 0.1 and A = 0.1, which corresponds to 
weak coupling in the adiabatic limit. In this regime the DMRG method gives less reliable 
results. The size of the polaron is underestimated, possibly due to finite-size effects in open 
boundary conditions. Note also that the DMRG method does not give symmetric results 
as it should, i.e. x(0 7^ x{~ 0- m ^ ms parameter regime our results have fully converged, 
as we can see from the perfect overlap of results of systems with two different sizes of the 
Hilbert space = 17, 18. Figure (^|c) plots correlations for uj = 0.1 and A = 0.435, which 
belongs to the strong coupling, small polaron regime. The DMRG produces superior results 
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in this regime, where our calculation at Nh = 21 has not fully converged to the large Nh 
limit. Our results are, nevertheless, in qualitative agreement with the DMRG. We conclude 
that both techniques give reliable results in the intermediate coupling regime, and that they 
complement each other in the weak and strong coupling regimes. 

A thorough investigation of correlation functions using ED and the variational Lanczos 
method was performed by Wellein and Fehske ||. Although we do not show a direct com- 
parison with their work, our results for correlation functions agree qualitatively with their 
calculations. 

While the strength of the DMRG calculation is exhibited in its ability to compute ground 
state properties of large systems, it is limited in its computation of excited states. In 
Figure (|^) we have shown how the nature of the polaron transforms from predominantly 
electronic character at k = to phononic around k = tt. We follow this transformation 
by computing the correlation function x f° r f° ur different values of k, shown in Figure ((7|). 
These parameters correspond to the weak coupling case in Figure (|5|). At k — 0, where 
the group velocity is zero, the deformation is limited to only a few lattice sites around the 
electron. It is always positive and exponentially decaying. At finite but small k = 7r/4, the 
deformation around the electron increases in amplitude and rings (oscillates in sign) as the 
polaron acquires a finite group velocity. At k = tt/2 the ringing is strongly enhanced. Note 
also that the spatial extent of the deformation increases in comparison to k — 0. The range 
of the deformation is maximum at k = tt, where it extends over the entire region shown 
in the figure. In keeping with the larger extent of the lattice deformation near k = tt, the 
ground state energy E(tt) converges more slowly with the size of the Hilbert space. 

We have also computed x f° r the strong-coupling case u = 0.8, A 2 = 3.2 (not shown). We 
find only weak k— dependence, which is a consequence of the crossover to the small polaron 
regime. The lattice deformation is localized predominantly on the electron site. 
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IV. WHAT IS THE NATURE OF THE FIRST EXCITED STATE? 



In this section, we focus on the question of whether an extra phonon excitation forms 
a bound state with the polaron, or instead remain as two widely separated entities. Using 
numerical and analytical approaches we will show that there exists a sharp phase transition 
between these two states. Although the ground state energy E is an analytic function of the 
parameters in the Hamiltonian, there are points at which the energy E% of the first excited 
state is nonanalytic. In previous work, Gogolin has found bound states of the polaron 
and additional phonon(s), but he does not obtain a phase transition between bound and 



unbound states because his approximations are limited to strong coupling \/uo > 1 p7 



A phase transition between a bound and unbound first excited state has been calculated 



for dimension d = 3 using a dynamical CPA approximation [28[ and dynamical mean field 
theory p9 . 



A. Numerical results 

We begin by computing the energy difference AE = E\ — Eq, where E\ and Eq are the 
first excited state and the ground state energies at k = (the two lowest bands in Fig. (^|) 
). In the case where the first excited state of a polaron can be described as a polaron 
ground state and an unbound extra phonon excitation, this energy difference should in the 
thermodynamic limit equal the phonon frequency, AE = u. In Fig. (|j) we plot the binding 
energy A = AE — uo for u = 0.5 as a function of the electron-phonon coupling A for various 
sizes of the variational space. We see two distinct regimes. Below A c ~ 0.95, A varies 
with the system size but remains positive (A > 0). Physically, for A < A c , the additional 
phonon excitation would prefer to be infinitely separated from the polaron, but is confined 
to a distance no greater than — 1 by the variational Hilbert space. As the system size 
increases, A slowly approaches zero from above as the "particle in a box" confinement energy 
decreases. In the other regime, A > A c , our data has clearly converged and A < 0. This is 



11 



the regime where the extra phonon excitation is absorbed by the polaron forming an excited 
bound polaron. Since the excited polaron forms an exponentially decaying bound state, the 
method already converges at Nh = 14. In the inset of Figure (§) we show the binding energy 
A in a larger interval of electron-phonon coupling A. Although the results cease to converge 
at larger A, we notice that the binding energy A reaches a minimum as a function of A. 
As we will demonstrate within the strong coupling approximation, the true binding energy 
approaches zero exponentially from below with increasing lambda. In Figure @ we show 
the phase diagram, valid for k = 0, separating the two regimes. The phase boundary, given 
by A = 0, was obtained numerically and using strong coupling perturbation theory in t to 
first and second order. Details of the latter calculation are given in the next subsection. 
The phase transition where A becomes negative at sufficiently large A is also seen in exact 



diagonalization calculations ||30|| . 



In Figure (|10|) we compute the distribution of the number of excited phonons in the 
vicinity of the electron 7(2 — j), Eq. (|7|), for the ground state 70 and the first excited state 71 
slightly below (A = 0.9), and above (A = 1.0) the transition for u = 0.5. The central peak 
of the correlation function 71 below the transition point is comparable in magnitude to 70 
(Figs. (]iTl|a,b)). The main difference between the two curves is the long range decay of 71 as 
a function of distance from the electron, onto which the central peak is superimposed. The 
extra phonon that is represented by this long-range tail extends throughout the whole system 
and is not bound to the polaron. See also the difference 71 — 70 in the inset of Fig. (|T0|b). 
The existence of an unbound, free phonon is confirmed by computing the difference of total 
phonon number A^i = X^7o,i(0- This difference should equal one below the transition 
point. Our numerical values give Nf h — Nq K ~ 1.02. We attribute the deviation from the 
exact result to finite-size effects. 

Correlation functions above the transition point (Figs. (|T0|c,d)) are physically different. 
First, phonon correlations in 71 decay exponentially, which also explains why the convergence 
in this region is excellent. Second, the size of the central peak in 71 is 3 times higher than 70 . 
(Note that to match scales in Figure (|i~0l d) we divided 71 by 3). The difference in total phonon 
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number gives 

N pn _ N pn _ 2 _ 33 _ We 

are thus facing a totally different physical picture: the 
excited state is composed of a polaron which contains several extra phonon excitations (in 
comparison to the ground state polaron) and the binding energy of the excited polaron is 
A < 0. The extra phonon excitations are located almost entirely on the electron site (see 
the inset of Fig. (|T0|d)). The value of 71 — 70 at j = is 2.16, which almost exhausts the 
phonon sum. 



B. Strong-coupling perturbation theory 

In this section, we calculate the ground and excited state energies of the polaron pertur- 
batively in the hopping t. For t = 0, the Hamiltonian Eq. (|l|) describes a harmonic oscillator 
with a shifted origin (due to the A force term) on the site with the electron, and unshifted 
oscillators on the other sites. The Hamiltonian in the new basis is given by the canonical 
transformation H = e s He~ s , where 

S = -gJ2 n j( a 3- a ])' ( 8 ) 
3 



and g = X/lu After some algebra (see also Ref. fll2|), the transformed Hamiltonian 



takes the following form: 



H = H + V 

H = uY^ a ) a j - ^g 2 n 3 
j 3 

V = -te~ 92 fctc i+ ie ff (°i- a i+i)e- a ^-^ +l) + h.c.) , (9) 



3 



where rij = c^Cj. The operator a] in Eq. @ creates a phonon excitation on site j relative to 
the shifted oscillator if there is an electron on site j, and relative to an unshifted oscillator on 
the other sites. The first term in Hq is the energy of the phonon excitations, and the second 
is the energy gained by the oscillator that is displaced by the force of the electron. In strong 
coupling perturbation theory, V in Eq. (g) is considered a perturbation. It represents the 
hopping of an electron, including possible creation and destruction of phonon excitations. 
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The lowest energy eigenstates of the unperturbed Hamiltonian H have no extra phonon 
excitations, and an energy Eq ^ = —X 2 /uj. They can be written \4>o(j)) = cj|(9), where 
| O) represents vacuum for electron and phonon degrees of freedom. This state represents 
a polaron localized on the site j. Evidently the ground state is N— fold degenerate, where 
N is the number of sites in the system. The perturbation lifts this degeneracy. The matrix 
elements V (i,j) = ((f> (i)\V\<f>o(j)) can be readily computed since the exponential factors in 
V are not effective in this case, 

V (t,j) = -te- 92 ; j = i±l. (10) 

This describes a translation-invariant tight-binding model in one dimension with nearest 
neighbor hopping. To first order in t, the ground state energy of a polaron with momentum 
k is 

E (k) = -\ 2 /uj -2te-' j2 cosk. (11) 

Excited states: The lowest energy states of H are the N degenerate states of energy 
—\ 2 /u, considered above. The next lowest energy sector contains the N 2 degenerate states 
of energy —\ 2 /cu + u, of the form = c^a\\0). The electron is on site j and the 

additional phonon excitation is on site /. We do degenerate perturbation theory to 0(t) 
in this excited sector. Using the total momentum k as a good quantum number, the 2d 
tight-binding problem in (j, I) becomes a Id tight-binding problem. The Id basis functions 
are |0i(j)) = \ipi(0,j)), where j is the distance between the phonon and the electron. The 
nonzero matrix elements Vi(i,j) = (<pi(i)\V\(f)i(j)) are 

14(0,0) = -2tg 2 e- g2 cos A; 
V 1 (0,±l) = -t(l-g 2 )e-° 2 
Vl (-l,l) = -tg 2 e- 92 e ik 

V 1 (i,j) = -te- 92 ; ./ / : 1:/../ / 0. (12) 

This is a Id tight-binding model that is translation-invariant except near the origin, where 
there is a second-neighbor hopping term and other modifications. The matrix Vi(i,j) is 
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Hermitian. Matrix elements in Eq. fll2|) define a secular equation \V(k) — E(k)\ = for the 
energies that can easily be solved numerically for a large system. For each total wavevector 
k, there are N independent solutions. The lowest energy first excited state at momentum k 
is 

E x (k) = -X 2 /uj + uj + E[ 1 \k), (13) 

where E^\k) is the lowest energy solution of the secular equation. 

The numerical and analytic solution of the secular equation reveals that there is a true 
phase transition (energy nonanalytic in the parameters) for the first excited state. This 
is perhaps surprising in light of the theorem that there can be no phase transition for the 



ground state ||24j| , and demonstrates the impossibility of extending the theorem to include 
excited states. We consider first a total momentum k = 0. For g > g x = 1, the lowest energy 
first excited state is found to be a bound state of a polaron and an additional phonon. The 
bound state is raman active. For g < g x , the first excited state is unbound, with an energy 
exactly u higher than the ground state. This energy is in fact an upper bound for the first 
excited state in any dimension and at any k, since one can construct a variational state 
with a zero momentum polaron and a momentum k phonon at infinite separation. The 
energy of the first excited state is nonanalytic (discontinuous first derivative) at g = g%. The 
bound state formation is unusual for several (related) reasons: (a) For g > g x , the binding 
energy is linear in (g — g x ) for small (g — g x ), rather than the (g — g x ) 2 dependence that is 
typical for Id bound state problems, (b) For g < gi, the phase shift is zero between the 
polaron and an additional unbound phonon of relative momentum q — > 0. The polaron and 
additional phonon pass through each other transparently, in contrast to the usual repulsive 
phase shift, (c) At g — g x , even for finite N (periodic boundary conditions), there are 
two exactly degenerate ground states, in contrast to the usual avoided crossing. These 
unexpected properties result in part from the fact that the central site in the tight-binding 
model (Eq. |T2|) becomes uncoupled from the rest of the lattice precisely at g = g x . The 
binding energy in Id can in fact be determined analytically from Eq. (0). One can show 
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that for g > gx, E[ 1 \k) = —te g ^(x 2 + l)/x, where 



3^-1 + ^-1)0^-1) . n 

x = ; k = 

2 

x = g 2 ; k = it . (14) 

The 0(t) strong-coupling analysis can be extended to higher dimensions and to nonzero 
total momentum k. The k = state always has an energy less than or equal to that at any 
other k. We find that for g < gi, the first excited state at any momentum k is unbound, 
with g\ — 1 in any dimension d > 1. The k = phase shift also vanishes for g < g\ 



in any dimension. For g > g 2 , the first excited state at any momentum is bound [32 
For g\ < g < g2i there is a phase transition on a surface in fc-space such that inside the 
surface (including k = 0) the first excited state is bound, and outside the surface (including 
k — (-7T, 7r, . . .) ) the first excited state is unbound. The location of the surface in fc-space 
is ^-dependent. The numerically obtained values are g 2 = 2.1 in 3d, g 2 = 1.66 in 2d, and 
92 — 9i — 1 m Id. (Id is special in that there is no intermediate phase. However, as g — > g± 
in Id, the k = it state is very weakly bound compared to k = 0; see Eq. (|i4|). The k = n 
state has zero amplitude to have the electron and the additional phonon on the same site.) 
The binding energy at k = is linear in (g — gi) for higher dimensions, as it is in Id. 
Strong coupling perturbation theory is carried to 0(t 2 ) in the appendix. 



V. CONCLUSIONS 

In summary, we have presented a variational approach for solving the Holstein model with 
dynamical, quantum phonons based on an exact diagonalization method. The variational 
space is defined on an infinite lattice. It is constructed by successive application of the off 
diagonal terms of the Hamiltonian starting from a single electron state. This leads to a 
systematically improvable variational basis that turns out to be efficient for calculating the 
ground state and low-lying excited state properties of the model. The method can compute 
properties in all parameter regimes, but it is at its best in the intermediate coupling regime, 
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where strong and weak coupling perturbation theories, and other variational methods have 
problems. 

The method allows the computation of energy bands and other physical properties at 
continuous wavevectors. For intermediate coupling strengths we are able to reach an accu- 
racy of 12 digits for the ground state energy at small k with as few as M = 88052 basis 
states and only a few seconds of CPU time on a workstation. Our results for energy bands 
presented in Fig. (H) compare well with other numerically more intensive methods and are 
more precise at both large and small k than some of them. Our energies are variational 
in the thermodynamic limit for any k. We believe that all results shown for energy bands 
computed with M = 3 x 10 6 states converge to at least 4 digits for arbitrary k. 

The accuracy of our method can be seen from the comparison of effective masses in 
Fig. (£|). While our results have converged to at least 4 digits in all parameter regimes 
presented in Fig. (||), GL and DMRG methods give less reliable results. While deviations 
of DMRG are insignificant, deviations of GL results near strong coupling seem to be sys- 
tematic. Results for the quasiparticle residue as a function of the wavevector in the weak 
coupling regime show a smooth crossover between the predominantly electronic character 
to a predominantly phononic character of the polaron. Our results agree with previous ED 
calculations. 

The correlation functions \ agree well with DMRG results in the intermediate coupling 
regime. In the weak coupling regime our method gives more reliable results. Close to 
the extreme strong coupling and adiabatic regime, our correlation functions have not fully 
converged as a function of for = 21, M = 6 x 10 6 . Our results in this regime are 
qualitatively similar to, but less accurate than those of the DMRG. Correlation functions 
X computed at different k provide detailed information on how the weak coupling polaron 
transforms as k increases from k = to k = 7r. 

Using numerical and strong-coupling approaches, we find a true phase transition (rather 
than a crossover) in the first excited state, where a polaron plus phonon system changes 
from unbound at weak coupling to bound at strong coupling. The first excited state does 
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not contribute to the optical conductivity, but rather is Raman active. 

There are a number of extensions to this work that we have not included in this paper. It 
is straightforward to consider anharmonic phonons. (In fact, the extreme double- well limit 
where only the lowest two states are retained is numerically less demanding than the linear 
case.) The AC conductivity and spectral function of a polaron can be calculated by the 
same methods. Properties of other Hamiltonians, including those with SSH-type couplings 
where phonons modify the hopping t can be calculated. Extensions that allow certain 
phonon excitations to be infinitely far from the electron are possible. Properties in higher 
spatial dimensions can be calculated. One can also calculate the properties of bipolarons, 
including Hubbard onsite and longer range interactions between electrons. And finally, one 
can calculate the coherent quantum dynamics of electron-phonon coupled systems driven 
far from equilibrium using similar methods |19| , |20[| . 

We would like to acknowledge valuable conversations with A. Alexandrov, A. Bishop, 
H. Fehske, E. Jeckelmann, L. C. Ku, K. Lindenberg, S. Marianer, F. Marsiglio, H. Roder, 
and S. White. J.B. would like to acknowledge the hospitality of Los Alamos National 
Laboratory where the major part of this work was performed. This work was supported by 
the US Department of Energy. 



APPENDIX: 



For simplicity we have limited our calculation of the energy corrections in second order 
strong-coupling perturbation theory to Id and zero momentum, k = 0. Following the work 
of Marsiglio the ground state correction to second order in the hopping t is given by 



E, 



(2) 



-It 1 



OJ 



E 



9 



2(m+n) 



+ 3E 



9 2n l 



n\ n 



(Al) 



=1 n\m\ n + m n _ i 

Calculation of the energy corrections of the excited state energy E\ involves degenerate 
perturbation theory, where matrix elements between degenerate states are computed 
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to second order in t. After a straightforward but tedious calculation we obtain for the 
non-zero matrix elements 



Vi(0,0) = -2tg 2 e 



+ 2t 



UJ 



g*W (l-^) 2 | ^ (l-^) 2 + 2 | 2 
-f^Lx n\m\ 1 — n — m ^ 2 n ' 1 — n 



^(0,1) = -t(l-0 2 )e 



- t 



UJ 



2(m+n) (1 _ _ ™) 

^ ' l-n-m ^ 



Ki(l,l)=t 2 



n,m=l 



n=2 



n! 1 — n 



+ 2 



a; 



E 



(7 2 ( m +") f g 2 {l-f) 2 i_ 

1 n!m! y 1 — n — m n + m / 



2n 2g i-2n+%(n + 2) 



n=2 



K 1 (_i > i) = _^ e -p' +f a2_!_ 



1 — n 

2n (1 _ » )2 



E 



n! 1 — n 



+ 1 



K(0,2) = -t 2 (l-^ 2 



_n=2 

29 2 q»i j 



n=l 



-V 



„2n i 

Mi,i + 2) = -* 2 — EV^ i>i 



w —j n! n 



V 1 (j,j) = -2t 2 



-2ff 2 



a; 



E 



g2(m+n) i g2n ^ 



+ 2E 



J>2 



(A2) 



(A3) 



(A4) 
(A5) 
(A6) 
(A7) 



(A8) 
(A9) 



\n,m=l Hlml U + m nTl Hl U _ 

V 1 (j,j + l) = -te~v 2 ; j>l. 
For k — 0, Vi (i, j) = Vi(j, z) = Vi(— i, — j). We numerically solve the secular equation 

(1 2 s ) 

\V — E\ =0. The lowest eigenvalue E\ of the secular equation gives us corrections to the 
excited state energy E\ to second order in the hopping t, 

Ex = -\ 2 /uj + uj + £f 1,2) . (A10) 
The binding energy to second order at k = is given by Eqs. ([O]), (|AT|), and ( |A10| ). 

A = E{^ 2) + 2te- g2 - E {2) . (All) 
The second order results are used to calculate the phase boundary in Fig. (§). 
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FIGURES 

FIG. 1. The small variational Hilbert space shown for the polaron is a subset of the = 3 
space. Basis states in the many-body Hilbert space are represented by dots, and nonzero 
off-diagonal matrix elements by lines. The x-coordinate of the dots is (aside from small displace- 
ments) the coordinate of the electron. Vertical bonds create phonons, and horizontal or nearly 
horizontal bonds are electron hops. State |1) is an electron on site and no phonons. State |2) is 
an electron and phonon, both on site 0. State |4) is an electron on site 1 and a phonon on site 0, 
which is reached from state |2) by hopping the electron to the right. State |5) is a translation of 
state 1 2). The Hamiltonian is sparse in this basis, with at most 4 bonds attached to a dot. The 
dots can also be thought of as Wannier orbitals in a one-body periodic tight-binding model. 

FIG. 2. The ground and excited state energy eigenvalues (those Ej < 0) are plotted as a 
function of k (in units of 7r) for A = oj = 1, Nh = 9, M = 1185. Excited states consist of the 
polaron with additional bound or unbound phonon excitations. The hopping parameter t = 1 is 
assumed here and throughout this work wherever t is not explicitly specified. 

FIG. 3. Polaron energy as a function of k (in units of it). Lines represent our results for 
two different sets of parameters obtained by calculating E{k) at 100 k-points. Circles, squares 
and diamonds are results obtained with Global-Local p5j , DMRG [15,17| and ED || calculations 
respectively. 

FIG. 4. The logarithm of the effective mass m* /mo as a function of X/u. Our results are 
plotted as full lines, and Global-Local results as dashed lines |l5|]. Open symbols, indicating the 
value of uj, are DMRG results IBT 



FIG. 5. The quasiparticle weight and in the inset the total number of phonons N^ h as a 
function of the wavevector k. Results by Wellein and Fehske || are represented by open circles. 
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FIG. 6. Lattice deformation x as a function of (i — j) for (a) u = 1.0 and A = 0.5, (b) u = 0.1 
and A = 0.1, and (c) uj = 0.1 and A = 0.435. Our results are represented by filled symbols and 
those of Jeckelmann and White |l7j by open circles. In case c) our results for Nh = 21 have not yet 
reached the large Nh limit. The star represents an extrapolation of our data (assuming exponential 
convergence) to Nh — ► oo, x(0) = 5.5 ± 0.1. 

FIG. 7. Lattice deformation x as a function of (i — j) for u = 0.8, A 2 = 0.4 and Nh = 18, for 
four different values of momentum k. The variational Hilbert space for Nh = 18 allows nonzero 
correlations to a maximum distance \i — j\ m ax = 17. Only distances up to 15 are plotted. 

FIG. 8. First excited state binding energy A = E\ — Eq — u as a function of A. Results are for 
uj = 0.5 and various Hilbert space sizes N^. Inset: binding energy over a wider range of A. 

FIG. 9. The phase diagram for the bound to unbound transition of the first excited state, 
obtained using the condition A(o>, A) = 0. The corresponding phase diagram for the ground state 
would be blank — there is no phase transition in the ground state, only a crossover. 

FIG. 10. The phonon number 7 as a function of the distance from the electron position (i — j) 
for the ground state (a) and the first excited state (b) both computed at A = 0.9 and the same 
in (c) and (d) for A = 1.0. All data are computed at phonon frequency u = 0.5 and = 18. 
Note that (d) is a plot of 71 /3. Insets in (b) and (d) represent differences 71 — 70 as a function of 
the distance (i — j). In (b), 71 — 70 drops to zero around \i — j\ = 15. This is a finite-size effect. 
Computing the same quantity with larger below the phase transition would lead to a larger 
extent of the correlation function indicating that the extra phonon excitation is not bound to the 
polaron. 
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TABLES 



TABLE I. 



\/uj 


Present 


ED N=16 


DMRG N=32 


Global-Local 


1 


-2.469684723933 


-2.46968477 


-2.46968 


-2.46931 


V2 


-2.998828186867 


-2.99882816 


-2.99883 


-2.99802 
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